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I. INTRODUCTION 

Modulational instability (MI) is a nonlinear phe¬ 
nomenon causing a plane wave (PW) or long pulse to 
break up into smaller, finite sub-structures. It can hap¬ 
pen in a plethora of different physical systems, ranging 
from deep water waves [I| to optical beams 0. In partic¬ 
ular, the mechanism of MI in optics, which may happen 
in both temporal Q and spatial domains [1, Ij, either 
independently or simultaneously, has been extensively 
studied due to its paramount implications on the dynam¬ 
ics of intense light beams or pulses traveling throughout 
nonlinear media. In brief, perturbations in the incoming 
optical field can be amplified upon propagation due to 
the nonlinear response of the optical medium 0 , leading 
to the reorganization of the energy and subsequent gen¬ 
eration of finite optical beams or filaments, randomly dis¬ 
tributed along the space-time profile of the parent wave. 

Regarding their effects, temporal MI of either contin¬ 
uous waves or very long pulses is responsible for the gen¬ 
eration of ultrashort pulses @ with temporal durations 
down to- or even below the single-cycle limit Q. Such 
Ml-driven short pulses are linked to the appearance of 
spectral sidebands which, in combination with other non¬ 
linear processes, can boost the generation of extremely 
broadband radiation. In contrast, spatial MI has been 
demonstrated to be, in general, an unwanted effect in 
the propagation of intense ultrashort pulses throughout 
bulk optical media. In fact, the Ml-triggered redistri¬ 
bution of the energy along the transverse profile of the 
beams gives rise to the appearance of intense hot spots 
that couldpotentially damage the optical material itself 
(see e.g. Q and references therein). Each intense spot 
contains an amount of power comparable to the criti¬ 
cal power Per for self-focusing Q- Whilst such struc¬ 
tures are doomed to undergo collapse due to self-focusing 
whenever their power P > P^r in a Kerr medium, they 
have been shown to form stable multidimensional solitary 
waves when traveling across optical media with compet¬ 
ing cubic and quintic (CQ) nonlinearities Q, yielding in 
some cases to the formation of liquid light states i.e., 


intense flat-top beams featuring intriguing surface ten¬ 
sion properties 0- Very remarkably, this new type of 
self-trapped beams have recently been observed experi¬ 
mentally in coherently-engineered atomic media [ll| . 

On the other hand, the recent measurement of Higher- 
Order Kerr (HOKE) responses in air and its con¬ 
stituents [l^ gave rise to an extensive, ongoing discussion 
on the physical origin and validity of such pec uliar con¬ 
tribution to the nonlinear polarization [13l-l20l|. Among 
other reasons, such a debate has been specially moti¬ 
vated by the deep implications of the HOKE response in 
the description of laser filamentation [2l| and novel light 
distributions [l^. In particular, theoretical media dis¬ 
playing such kind of competing HOKE terms have been 
shown to allow for a new branch of highly concentrated 
nonlinear solitary waves, called ultrasolitons [^ . which 
might coexist with the usual solitons appearing in CQ- 
like nonlinear models [^ . 

In this paper, we introduce a complete analytical and 
numerical study of the MI process in a system governed 
by a canonical nonlinear Schrddinger equation (NLSE) 
involving a local, arbitrary nonlinear response F to the 
applied field. In particular, in Sec. H we will revisit 
the theory introduced in Q to account for the recently 
proposed HOKE nonlinearities, deriving very simple ana¬ 
lytical criteria for the identification of multiple regimes of 
stability and instability of PW solutions in such systems. 
To illustrate our results, in Sec. HI we will show that the 
new model explains the phenomenology reported in the 
literature for different kinds of nonlinear responses. In 
addition, we will discuss the parametric regime of HOKE 
responses proposed in , which allows for the existence 
of ultrasolitons. In this case, we will show that the theory 
predicts several alternating stability-instability windows 
whose existence is also demonstrated by means of nu¬ 
merical simulations in Sec. IV. Thus, our results define 
a completely new MI landscape as compared to the sce¬ 
nario that has recently been described [25l | in the presence 
of HOKE nonlinearities out of the multistability regime. 
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II. GENERAL THEORY 

A. The canonical scalar NLSE and its plane wave 
solutions 

We consider a system described by the wave function 
dt (r, rj) evolving along the rj direction. Such a system can 
be either a nonlinear optical medium in which dt(r, 77 = z) 
would represent the scalar electric field envelope of an 
optical wave propagating along the 2; direction , or an 
ultracold atomic gas in which 4'(r,77 = t) describes the 
order parameter of the corresponding macroscopic col¬ 
lective quantum state evolving in time t [2^. In such a 
systems, the evolution of \E'(r, rj) in a (n-b l)-dimensional 
space of points (r, ry) is governed by the canonical non¬ 
linear Schrodinger equation 

2 ^ + = 

where accounts for diffraction (dispersion) 

in the spatial (temporal) domain and F is an arbitrary, 
real-valued and continuous function of |'I'|. For the sake 
of simplicity, we will restrict our discussion to Eq. o, 
although our approach can be easily generalized in the 
presence of terms involving higher order derivatives, such 
as those that should be introduced to describe the prop¬ 
agation of ultrashort pulses [2^. 

We assume that all the quantities in Eq. o, including 
the coordinates r, have been suitably rescaled in such a 
way that they are dimensionless. All the relevant param¬ 
eters will then appear in the explicit form of the function 
F. 

It is well-known that Eq. © admits PW solutions of 
the form 

vk = A exp(-i/ir7 -b iipo), (2) 

where A > 0 , fi and (po are real constants describing 
the amplitude, propagation constant (chemical potential 
in case of matter-waves) and global phase, respectively. 
Eq. © implies the relation 


B. Modulational instability and perturbation 
growth rate 

Let us now study the stability of the PW solutions of 
Eq. ([ 1 } under the influence of small perturbations by 
generalizing the perturbative method of Ref. 0 . Simi¬ 
lar analyses have also been carried out in Refs. [27l - [^ . 
We will then look for a solution of the form 'k (r, 77) = 
(A -b /(r, 77)) exp(—77x77 -b 7(730)1 where f{v,r]) is an arbi¬ 
trary, complex-valued perturbation. An additional in¬ 
teresting possibility, which lies beyond the scope of the 
present paper, would correspond to the consideration of 
localized perturbations on finite backgrounds [s^ . [3^ . 

The requirement that the pure PW solution is stable 
under the action of small perturbations corresponds to 
the condition that |/(r, 77)] can be kept |/(r,77)| <g; A for 
all values of r and 77. In this case, at the first order in /, 

= |A -b ff = A^ -b 2Aa + a^+(3^ - A^{1 + 2a/A), 

where a and /3 are real functions describing the real 
and imaginary parts of /, such that f{r,r]) = a{r,ri) + 
7 / 3 (r, 77 ), and |a|,|/ 3 | <C A. Therefore | 4 '| ~ A + a, and 
Eq. © becomes 

% + + <^AF\A) = 0 , ( 4 ) 

where F'{A) = ^^(A) and we have taken into account 
Eq. ©. By separating the real and imaginary parts, we 
get a set of two coupled equations, 

-% + \y‘^a + AF'{A)a = {). ( 5 ) 

This system of equations can be solved in the reciprocal 
space by expanding Q;(r, 77) = Q!(k, 77) exp(7k • r)(i"k 

and / 3 {r,r]) = ,d(k, 77) exp(7k • r)d"k, where (a, / 3 ) 

stand for the Eourier transforms of (a, P) and n is the 
dimension of the space of the vectors r. We then obtain 
P = being k = |k|, together with the following 

equation 


fi = -F{A), ( 3 ) 

which reflects the nonlinear behavior of the system since 
the phase of the wave is self-modulated by its intensity Q . 
In addition, due to the U(I) symmetry of Eq. ([T]), its dif¬ 
ferent solutions are invariant under global phase transfor¬ 
mations and so the arbitrary phase factor ipQ will play no 
role in the following discussion. 

In the next sections, we will analyze the dynamics of 
PWs propagating through optical media governed by Eq. 
© in the presence of noise. In particular, we will derive 
a fairly simple analytical rule for the prediction of the 
onset of MI. Such a rule is completely general, and can 
be applied to a plethora of nonlinear systems described 
by different nonlinear responses F. 


^ -{-uj^{k)a = 0 , ( 6 ) 

which is formally similar to the equation governing the 
dynamics of a classical harmonic oscillator whose angular 
frequency uj{k) satisfies the following dispersion relation 
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AF'(A) 


( 7 ) 


Whenever uj{k) becomes pure imaginary for a certain 
value of k, the perturbations described by (a, j 3 ) experi¬ 
ence an exponential amplification. In this case, the cor¬ 
responding PW will be unstable and the perturbative 
approximation will be no longer fulfilled. On the other 
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hand, if uj{k) is real for all wavevectors fc, the values of 
a and /3 will be constrained by the initial conditions. In 
other words, the perturbations f{r,ri) can be kept small 
at all points of the space (r, 77 ) if and only if uj{k) is a 
real number for all values oi k. As a consequence, the 
stability condition can be expressed in the following very 
simple form 


F\A) < 0. 


( 8 ) 


Accordingly, whenever F'{A) > 0 the corresponding 
PW solution will not be stable and will eventually break 
up into small localized substructures or filaments 
In the latter situation, it is very useful to compute the 
wavevector fcmax related with the fastest growing Fourier 
component by maximizing the value of the exponential 

1/2 


growth factor r(A:) = iuj{k) = AF'{A) - ^ 
obtain 


. We 


F(|vk|) = (11) 

and the constant /2 can be chosen to be ±1 by suitably 
rescaling 'I'. The stability condition Eq. ([5]) can then be 
written as /2 < 0 , corresponding to the self-defocusing 
Kerr nonlinearity. In the opposite case, /2 > 0, the PWs 
will undergo MI for any value of the amplitude A. The 
fastest growing perturbations leading to the break-up of 
the plane wavefront in multiple filaments will correspond 
to a wavevector 


fcmax = \/2A, (12) 

and the corresponding value of the maximum perturba¬ 
tion growth rate will accordingly be 


fcmax= [AE'(A)]'/^ (9) 


r — 4 ^ 

max — • 


(13) 


This will be our theoretical prediction for the location 
of the peaks of the MI sidebands, as they are called in 
optics. The corresponding prediction for the maximum 
exponential growth rate, called maximum MI gain in the 
context of nonlinear optics, will then be 


As expected, these results agree with those obtained 
in the original paper by Bespalov and Talanov Q. 


B. Cubic-Quintic nonlinearity 


Pmax — r(^max) — 


AF'{A) 


( 10 ) 


For F'{A) > 0, and in the presence of noise, the per¬ 
turbed PW will break into filaments having a character¬ 
istic size given by A = Tr/fcmax = tt/ \/2Vraax 0- Notice 
that the expression for A is only valid at the very first 
stages of the wave destabilization. 


CQ materials have been thoroughly studied in the lit¬ 
erature from the theoretical point of view due to 
the interesting properties they display as a result of their 
nonlinear response 00 - Remarkably, the first real¬ 
ization of a canonical CQ nonlinearity has recently been 
accomplished in coherent atomic optical media [llj. 

By suitably rescaling t[/ and the space-time coordi¬ 
nates, light propagating throughout CQ systems can be 
described by Eq. o involving a nonlinear term 


III. EXAMPLES OF APPLICATION. HOKE 
NONLINEARITIES. 

In this section, we will apply the general formalism 
derived above to some canonical nonlinear systems like 
the Kerr and Cubic-Quintic (CQ) nonlinear media Q. 
In particular, we will show that the general expression of 
Eq. ([ 8 |) reproduces the known results on the MI dynamics 
for such media. In addition, we will discuss the case of 
a generalized HOKE nonlinearity, and we will find the 
appearance of different instability windows, depending 
on the relative strengths of the higher order nonlinear 
terms. 


=/2(I^P - I'l^l") (14) 

where = ±1. 

Taking into account that we have chosen A > 0, the 
stability condition Eq. ([5]) can then be written as /2(1 — 
2A2) < 0, i.e. A > 1/^/2 for /2 = -bl > 0, or A < 1/^2 
for /2 = —1 < 0. 

Again, the opposite case, /2(1 — 2A^) > 0, i.e. A < 
l/-\/2 for /2 = -1-1, or A > l/\/2 for /2 = —1, will yield 
to MI of the incoming PWs, and the fastest growing per¬ 
turbations will then feature a wavevector 


A. Cubic nonlinearity 

The cubic nonlinearity (called Kerr nonlinearity in op¬ 
tics) naturally appears both in the mean-field description 
of Bose-Einstein condensates and in the modeling of a 
large class of optical systems. In this case. 


^max — y2Av'|l-2A2|, (15) 

being the maximum perturbation growth rate 

rmax= |A2-2A4|. (16) 

These results agree with those obtained in Ref. 0 . 
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FIG. 1. (Color online) Multistability regions for HOKE mod¬ 
els in the (/ej/s) parameter space. The inner (pink) area is 
described by Eq. (|20l) and corresponds to the existence of soli- 
tonic multistability (see Ref. [23]). The wider (purple) area 
corresponds to the existence of two branches of stable PWs, 
as described by Eq. IISI). 


As we have chosen A > 0, the stability condition Eq. 
dSj) can be written as 1 — 2V + < 0, where 

V = A^. We can then have two different scenarios, de¬ 
pending on the roots of the algebraic equation 

1-2V + - 4faV^ = 0, (18) 

namely: i) if Eq. (ITS)) has only one positive real root 
Vi, the stability condition reads A > \fV\, while for 
A < '/V\ the PWs will undergo MI. ii) If Eq. (IT^ has 
three positive real roots Vi < V 2 < Ps; we find two dif¬ 
ferent stability windows: \fV\ < A < sfV^ and A > 
whereas for A < \fV[ and for sfV^ < A < y/Vs the 
corresponding PW will be modulationally unstable. We 
can obtain an analytical condition on the parameters /g 
and /g for Eq. m to have three real roots, which will 
delimitate the four stability and instability regions dis¬ 
cussed above. This can be accomplished by calculating 
the {fe, fa) pairs for which the discriminant associated 
with the polynomial expression Eq. (ITS)) turns ont to be 
negative, giving 


27/1 - 9/1 - lOS/g/g + 32/g + 108/1 < 0. (19) 

C. Higher-Order Kerr nonlinearity 


HOKE nonlinearities have recently been introduced 
to describe the experimental observation of a saturation 
(and subsequent sign inversion) of the nonlinear correc¬ 
tion to the refractive index at high optical intensities in 
gases M- Among other implications, they have been 
theoretically argued to provide a new mechanism for the 
stabilization of optical filaments [Slj, js^ without the need 
for plasma-related effects [2l| . Such a nonlinear response 
can also give rise to the existence of new localized light 
structures, called fermionic light states MM, liquid 


light states [ilia, Si, [231, or ultrasolitons |23j| . By suit¬ 
ably rescaling if and the space-time coordinates, we will 
model a general nonlinear HOKE response as 



= (17) 

where, to be concrete, we will consider f 2 q > 0 for all q 
and /4 = /2 = 1 like in the aforementioned CQ case. In 
the following discussion, motivated by the measurements 
in air and oxygen [liSil, we will also assume that n = 4, 
so that f 2 q = 0 for g > 5. 

The dimensionless parameters entering the NLSE are 
related to the refractive index n = no + An = no + 
'Yliq=i'^ 2 ql‘^, where no is the linear refractive index, by 
the relations An = (n|/|n 4 |)E', ng = {n1/n2)fe and 
n-a = (ii 4 /ii 2 )/ 8 , as shown in Ref. [l^. The only free 
parameters included in Eq. [I]are then /g and fa, which 
will be assumed to be positive taking into account the 
results of Ref. [l^ for the optical response of common 
gases. 


FIG. 2. (Golor online) Representation of Vmax vs. A for the 
HOKE model with /g = 2.8 and fa = 3.9. Whenever posi¬ 
tive, Fmax gives the exponential growth factor of the fastest 
growing filaments. 

As shown in Fig. [1] this condition (purple area), allow¬ 
ing for the existence of two different parametric regions 
for stable PWs, is less stringent than the condition for 
the existence of multistable flat-topped soliton solutions 
reported in Ref. (pink area), namely 

18225/1-5400/1-77760/g/8-h 204802/8-h 93312/1 < 0. 

( 20 ) 

Hence, whenever two different soliton branches ex¬ 
ist, thus implying soliton multistability as discussed in 
Ref. [ 2 ^ , there will also be two different regions of stabil¬ 
ity for the PWs. However, the existence of two stability 
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FIG. 3. (Color online) Representation of Tmax vs. A for the 
HOKE model with fe = 0.3 and fs = 0.02. Whenever posi¬ 
tive, Fmax gives the exponential growth factor of the fastest 
growing filaments. Inset: close-up of the low-amplitude re¬ 
gion. 


regions for the PWs does not necessarily imply the emer¬ 
gence of a second, ultrasolitonic branch for the localized 
solutions. 



FIG. 4. (Color online) Evolution of the peak amplitude of 
a perturbed PW with A=0.55 in a HOKE medium {fe = 
0.3, fs — 0.02). The results of the numerical simulation are 
represented by the black line, whereas the red line corresponds 
to a fit of the evolution to Eq. (1231) . The coefficients of the 
analytical fit are displayed in Table 2. 

Following the same procedure described above, for any 
value of A belonging to one of the instability windows 
discussed above, we can compute the wavevector corre¬ 
sponding to the fastest growing perturbations as 


/Cmax = . (21) 

being the corresponding value of the maximum pertur¬ 
bation growth rate 

n 

r„,ax = (22) 

9=1 



Simulation 

Theory 

A 

a (xlO^®) 

b 

c 


■^sim 

r 

max 

A 

0.30 

2.784 

0.090 

0.304 

0.9998 

8.15±0.63 

0.079 

7.91 

0.60 

29.704 

0.210 

0.607 

0.9993 

4.93±0.27 

0.231 

4.62 


TABLE I. Results of the numerical simulations for the evolu¬ 
tion of perturbed PWs with different amplitudes in the HOKE 
model with fe = 2.8, fs = 3.9. The first column gives the 
value of the initial amplitude A. The next columns give the 
values of a, b and c obtained from the fit of the evolution 
of the PW peak amplitude to an exponential growth, as de¬ 
scribed in Fig. m gives the correlation of the fit. Asim 
is the numerical estimate for the initial width of the fastest 
emerging filaments. The last two columns reproduce the the¬ 
oretical predictions of Tmax and A = 'k f \/2Vmax for compari¬ 
son. Notice that the parameter h is the numerical counterpart 
of F 

max • 


Let us illustrate these results with two concrete exam¬ 
ples: 

i) /e = 2.8 and /g = 3.9, corresponding to the central 

values of the n 2 q obtained in the experiment of Ref. [l^ 
to describe the propagation of ultrashort laser pulses in 
oxygen, namely n 2 = 1.6 x 10 “^®cto^/FF, n 4 = —5.2 x 
lt)-^^crrA/W'^, ne = 4.8 x /W^, ns = -2.1 x 

10“®®cm®/IF^. In this case, Eq. (fTOl) is not satisfied, so 
that Eq. (IT51) has only a single real root, Vi = 0.526, 
corresponding to the amplitude Ai = ^/Vl = 0.725. The 
dependence of F^ax on the PW amplitude A, as given 
by Eq. is plotted in Eig. [2] Notice that above 

the threshold A > Ai, Fmax < 0 indicating that the 
corresponding PWs are linearly stable. 

ii) /e = 0.3, fs = 0.02, corresponding to the param¬ 
eters introduced in Ref. to discuss a HOKE model 
predicting soliton multistability (Eq. (OUl) is fulfilled in 
this case). Moreover, with these parameters Eq. m is 
also satisfied, and so we get three real roots of Eq. m, 
namely Vi = 0.716, V 2 = 2.060 and V 3 = 8.474, whose 
square roots give the amplitudes Ai = 0.846, A 2 = 1.435 
and As = 2.911, respectively. The corresponding depen¬ 
dence on A of the values of Fmax, as computed through 
Eq. (l22l) . is plotted in Eig. 0 As we can clearly see in the 
figure, there are two stability windows, where Fmax < 0 , 
for Ai < A < A 2 (see the inset of Fig. [31) and for A > A 3 . 
Conversely, there are two instability regions spanning the 
parametric ranges A < Ai and A 2 < A < A 3 , respec¬ 
tively. 

Interestingly, it has been recently reported that the 
HOKE response is usually non-instantaneous and it may 
also show a further effective dep endence on both the in¬ 
tensity and pulse duration [T^ j^. Thus it is worth 
investigating different ranges of the {fejs) parameters, 
besides the values reported in [l^. Our Eig. 1 can be 
used to infer how the effective parameters fe and fs vary 
without significantly modifying the qualitative behavior 
of the propagation of light in the media. Indeed, our 
analysis indicates that similar modulational instability 
landscapes are expected for values of the parameters fe 
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and fs lying within the same region of Fig. 1. On the 
other hand, this may also open the possibility of the ex¬ 
istence of multistability and “ultrasoliton” states even 
in common optical media, provided that the effective 
(/sj/s) values enter in the multistability region for cer¬ 
tain suitable intensity levels. At this respect, ab-initio 
and perturbative calculations for the hydrogen atom have 
revealed that an effective HOKE response can only be tol¬ 
erated when all the higher-order corrections correspond 
to about 10% of the cubic response [s^. Under those con¬ 
straints, we have checked that a sensitive HOKE response 
could be observed in a hypothetical medium featuring 
/e = 0.3, fs = 0.02 for intensities above 18 TW/cw?. 

In the next section we will illustrate the theoretical 
results reported above by showing some numerical ex¬ 
amples of the evolution of PWs in different nonlinear 
systems. 


IV. NUMERICAL SIMULATIONS 

Henceforth we will discuss the results of our numerical 
simulations for the evolution of the systems described in 
the previous section when different initial PWs are per¬ 
turbed with low-amplitude white noise. We employed 
a standard split-step beam propagation method to solve 
Eq. ([T|) , being the step involving the linear part of the 

operator (diffraction) treated using a hnite-differences 
scheme, the so-called Crank-Nicolson algorithm [sgI - I^ . 
together with Neumann boundary conditions ^ = 0 at 
the borders. As this method makes efficient use of the 
whole computational domain for the physical results, we 
can mimic the evolution of PW-like fields featuring hnite 
energy with an affordable computational effort. There 
are other approaches to simulate PW-like beams such as 
the consideration of large, but finite, flat-top beams re¬ 
stricting the observation plane to the central part of the 
domain, or the usage of periodic boundary conditions. 
These techniques, however, often introduce non-physical 
interactions in the propagation, which may eventually al¬ 
ter the results. Thus, in our case we can safely consider 
wide enough PWs and study the evolution of central ar¬ 
eas of width much larger than the characteristic radii of 
the filaments born after the PWs destabilisation. 

By monitoring the exponential growth of the maxi¬ 
mum amplitude of the unstable fields, we can numeri¬ 
cally obtain a direct estimation of the largest perturba¬ 
tion growth rates. To do so, we will fit the evolution of 
the maximum amplitude |'I'(77)| during the early stages 
of the PW destabilization to an exponential of the form 


|^|=ae'”' + c, (23) 

where the theoretical values for b and a -I- c correspond to 
Tmax and A = |'I'(0)|, respectively. An example of this 
exponential behaviour of |^'(?7)| is shown in Eig. 31 where 
we represent the early stages of evolution of a PW with 
initial amplitude A=0.55 in a HOKE nonlinear medium. 



FIG. 5. (Color online) Simulation of the evolution of a PW 
of initial amplitude A = 0.3 (perturbed by white noise) in 
a system displaying HOKE nonlinearity with fs = 0.3 and 
fs = 0.02. The different panels show the two-dimensional 
amplitude distribution |^(x,t/)| of the incoming field at dif¬ 
ferent stages of its evolution (a) rj = 80, (b) rj = 115, (c) 
rj = 125 and (d) rj — 150. The spatial region displayed in the 
figures corresponds to x,y £ [100,100]. 


together with a fit of the simulation results to Eq. (1231) . 
As it can be inferred from the graph, the agreement be¬ 
tween simulation and theory is remarkable. 

As the MI processes for the Kerr and CQ nonlinear 
models have already been extensively studied in the lit¬ 
erature, we will concentrate the following discussion to 
the HOKE model for two different parametric regimes, 
namely (/e = 2.8, fs = 3.9) and (/e = 0.3, fs = 0.02). 


A. HOKE Model with fs = 2.8 and fs = 3.9. 

Let us consider a system displaying a HOKE nonlin¬ 
earity of the form described in Eq. (33 with /e = 2.8 
and fs = 3.9. In this case, as we have discussed above 
and as shown in Fig. [2l the theory predicts the occur¬ 
rence of MI for initial amplitudes A < Ai = 0.725, and 
stability for A > Ai. The physical mechanism behind 
the stabilization of the PWs with amplitudes above the 
estimated threshold Ai is the following: for small or mod¬ 
erate values of A < Ai, the positive terms entering the 
effective nonlinear refractive index correction F (namely 
those proportional to /2 and /e) dominate, and hence F 
increases with increasing A. This regime has been shown 
to support the existence of localized solitary waves sat¬ 
isfying an equation of state formally similar to that gov¬ 
erning the dynamics of a degenerate Fermi gas . For 
values of A > Ai, however, the negative terms entering 
F (namely those proportional to f^ and fs) dominate, 
and the function F decreases with increasing A. Thus, 
in analogy with both Kerr and CQ systems, small per- 
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turbations on a PW background cannot be developed 
into hot filaments when the effective nonlinearity is self- 
defocusing, since light will tend to ffow away from the 
regions of high intensity. 



0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 


FIG. 6 . (Color online) Simulation of the evolution of a PW 
of initial amplitude A = 2.55 (perturbed by white noise) in 
a system displaying HOKE nonlinearity with fe = 0.3 and 
/s = 0.02. The different panels show the two-dimensional 
amplitude distribution |^(a;,i/)| of the incoming field at dif¬ 
ferent stages of its evolution (a) 77 = 0.3, (b) r] = 0.5, (c) 
ry = 1 and (d) rj = 5. The spatial region displayed in the 
figures corresponds to x,y & [10,10], since the characteristic 
size of the emerging filaments (see panel (a)) is much smaller 
than that of the filaments displayed in Fig. [5] 

In this framework, the results of our simulations for the 
evolution of two perturbed PWs with amplitudes A = 0.3 
and A = 0.6, i.e. well within the instability window, are 
summarised in Table 1. We have corroborated that the 
maximum amplitude evolves according to the exponential 
law given by Eq. (1^51) in a very good approximation. In 
particular, the numerical growth rates b calculated from 
the fits to Eq. agree with the theoretical predictions 
^max obtained from Eq. (|22|) . Furthermore, we have also 
compared the analytical estimations of the characteristic 
size of the arising filaments A with the numerical results. 
To do so, we have monitored the size of the amplitude 
modulations Agim by subtracting the initial PW back¬ 
ground at the very early stages of the PW destabilization. 
The results of this analysis are also included in Table 1, 
where the numerical value of Agim has been obtained by 
averaging over an ensemble of several filaments. 


B. HOKE Model with fe = 0.3 and /s = 0.02 

In this section we will analyze an example of HOKE 
model involving two instability windows. We choose 
/e = 0.3, /§ = 0.02, which lies within the multistabil¬ 
ity region as shown in Fig. |T] above. In particular, we 


will study the dynamics of two PWs with initial ampli¬ 
tudes A = 0.3 and A = 2.55, belonging to the first and 
second instability regions introduced above, respectively 
(see section III-C and Fig. |3|). Figs. [5115] summarize 
the main results of the simulations for the evolution of 
these two PWs. In Fig. (U we observe that the initially- 
perturbed PW featuring A = 0.3 undergoes a redistri¬ 
bution of its energy as it evolves (panels (a)-(b)), giving 
rise to the formation of soliton-like localized structures 
(panel (c)). These solitary waves can then interact with 
each other in a complex way, provided that their phases 
are not correlated. Such dynamical behavior in nonlinear 
media has been first demonstrated in [H, . 

After a certain evolution period in which all emerg¬ 
ing nonlinear beams exchange energy among them, some 
of the structures stabilize and form perturbed 2D soli- 
tons (panel (d)). Interestingly, all remnant solitons fea¬ 
ture Gaussian-like shapes resembling those of the solitary 
waves found in Ref. |22| for moderate amplitudes, i.e., far 
from the strong self-defocusing regime where fiat-topped 
beams do exist [ 13 . In particular, their amplitudes are 
slightly above the limiting value Af^^ = 1.09, which cor¬ 
responds to the asymptotic PW solution of the ordinary, 
CQ-like soliton branch (2^ . This asymptotic behavior is 
reproduced in Fig. 0 As we can appreciate in the figure, 
after a certain transient period the maximum amplitude 
of the field stabilizes, which turns out to coincide with 
the emergence of stable solitary waves in the system as a 
consequence of the Ml-driven break up process. 



FIG. 7. (Golor online) Evolution of the maximum amplitude 
of the PW fields with initial amplitudes A = 0.3, A — 0.55, 
belonging to the first instability window of the multistable 
HOKE model involving /e = 0.3 and /g = 0.02. Solid lines 
correspond to the numerical simulations, whilst the dashed 
line represents the limiting PW solution of the localized soli¬ 
ton branch of the system. 

On the other hand. Fig. |6]shows the dynamics of a per¬ 
turbed PW featuring A = 2.55, i. e. well within the sec¬ 
ond instability window displayed in Fig. O Again, we ob¬ 
serve a redistribution of the energy in the PW front with 
the subsequent formation of localized structures (panels 
(a)-(b)). However, in contrast to the behavior shown 
in Fig. O the perturbed solitons emerging during the 
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b 
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A" 

^sim 

r rnax 

A 

0.30 

5.388 

0.085 

0.304 

0.9998 

8.22±0.28 

0.075 

8.14 

0.55 

9.766 

0.149 

0.557 

0.9998 

6.02±0.23 

0.144 

5.86 

0.75 

12.517 

0.080 

0.760 

0.9974 

8.20±0.81 

0.082 

7.76 

1.60 

0.556 

1.230 

1.602 

0.9997 

2.26±0.24 

1.116 

2.10 

2.55 

9.713 

26.197 

2.553 

0.9999 

0.53±0.05 

26.361 

0.43 


TABLE II. Results of the numerical simulations for the evolu¬ 
tion of perturbed PWs with different amplitudes in the HOKE 
model with /e = 0.3, fs = 0.02. All parameters displayed are 
the same as those included in Table 1. 



FIG. 8. (Color online) Same as in Fig. [7] with initial ampli¬ 
tudes of the fields A — 1.6 and A = 2.55, belonging to the 
second instability window of the multistable HOKE model 
involving /e = 0.3 and fs = 0.02 


dynamics merge together to form large-scale structures 
with a nearly-homogeneous amplitude, as it can be in¬ 
ferred from the fields displayed in panels (c)-(d). This 
suggests that, for A = 2.55, we can reach the threshold 
of existence of flat-topped states Hi 113) nil , so that the 
system evolves asymptotically towards the formation of a 
different PW-like structure, whose particular amplitude 
can be estimated analytically to be = 3.21 [^ . 

Again, the stabilization of the system at values slightly 
above is corroborated by monitoring the dynamical 
evolution of the peak amplitude, as displayed in Fig. [S] 

In addition, it is worth noticing that the overall dy¬ 
namics displayed in Fig. [5] turns out to be much faster 
than that represented in Fig. [51 The reason is that Tmax 
peaks close to A = 2.55 according to the theory (see Fig. 


El), and its absolute value is more than two orders of mag¬ 
nitude higher than that corresponding to A = 0.3, thus 
justifying the much shorter destabilization scale observed 
in Fig. [51 In order to reinforce this assertion, we have car¬ 
ried out the same numerical simulations described above 
with PWs featuring amplitudes A = 0.55 and A = 1.6, 
whose results are also depicted in Fig. 0 and Fig. [51 
respectively. We observe that, as expected, the destabi¬ 
lization occurs earlier since, in this range, Tmax increases 
for larger amplitudes following Eq. (1^^ . In light of the 
same procedure described in Section IV-A above, we have 
also numerically computed Tmax by fitting the exponen¬ 
tial growth of the P W amplitude at the early stages of MI 
to Eq. (1^51) . The estimations of Fmaa,, as well as the char¬ 
acteristic size of the fastest growing filaments Agim, are 
summarized in Table 2 for different initial values of the 
PW amplitude. One can appreciate that the agreement 
between the analytical and numerical results is reason¬ 
able, even though there are small discrepancies inherent 
to the the fact that the first-order perturbation theory is 
only applicable when the amplitudes of the perturbations 
are very small. 

V. CONCLUSIONS. 

In this paper we have put forward a full analysis of the 
nonlinear process of modulational instability in systems 
described by nonlinear Schrddinger equations with arbi¬ 
trary instantaneous nonlinear responses. In particular, 
the proposed theoretical approach allows for a complete 
description of the multiple regimes of stability and in¬ 
stability of plane waves in systems involving competing 
higher-order Kerr nonlinearities. All our analytical pre¬ 
dictions for the stability domains, perturbation growth 
rates and characteristic filament sizes have been con¬ 
firmed by direct numerical simulations of the evolution 
of unstable plane waves. This intriguing phenomenol¬ 
ogy could potentially be observed in a coherent atomic 
medium with a properly tailored nonlinear response. 
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